3 Variant Calling

3 Variant Calling

contact: ludwig@univ-brest.fr

Principles

Once the Reads are aligned against a Reference Genome the Variant Calling step will look for differences between the reads and the Reference Genome at each genomic position.

variantcalling.png

3.1. Tools

Several Algorithms for Germline Variants

Algorithm True Positive False Positive
HaplotypeCaller ⭐ 21,707 367
FreeBayes 23,857 18,256
SAMtools Mpileup 25,081 2,129
SNPSVM 17,920 65

Algorithms for Somatic Variants

  • MuTect2
  • Strelka
  • EBCall
  • Seurat
  • VarScan2

3.2. Generating a gVCF with Haplotype Caller (GATK4)

gatk HaplotypeCaller \
	-I child.bqsr.bam \
	-R ref.fa \
	-L /DATA/bioinfo/chr12M.capture.bed \
	-ERC GVCF \
	-O child.gvcf.gz

bash
Parameters explained

-I input bam file
-R reference genome
-L positions covered by the Capture Kit
-ERC GVCF outputs intermediate GVCF
-O output file name

VCF vs. gVCF file formats
  • VCF files provide information only for sites where at least one sample has a variant allele
  • gVCF files (as produced by HaplotypeCaller) provide a genotype for each position and not just variant positions. This is useful to differentiate missing calls from position homozygous to the reference genome.
Regions: BED Format

Sometimes we want to work on a subset of variants, defined by regions. Those regions can for example match a gene/set of genes or the positions taken into account by the capture kit.

The prevaling format to define regions is the BED (Browser Extensible Data) file format. It a tabulated file, with the following columns:

# Name Description
1 chrom Chromosome (e.g. chr3, chrY, chr2_random) or scaffold (e.g. scaffold10671) name
2 chromStart Start coordinate on the chromosome or scaffold for the sequence considered. In base
3 chromEnd End coordinate on the chromosome or scaffold for the sequence considered. In base
  • base means that the first nulceotide of the chromosome is 0
  • base means that the first nulceotide of the chromosome is 1 (can also be view as base excluded : [start;end[)

VCF File are numbered in base, so if you want the chromosome 2 on position 11,12,13,14 et 15 you don't write chr2 11 15 but chr2 10 15 and for just the position 17, you write chr 16 17

3.3. Generate Data for mother and father

At this time only the files for the child are available. Copy the files for the mother and father (and files you may have missed for the child).

cd
/DATA/bioinfo/getmissing.sh

bash

3.4. Combining multiple gVCF files

  • gVCF files produced by HaplotypeCaller contain a single sample
  • The actual genotyping (next step) will need to look simultaneously at the genotype of each sample
  • If a large (200+) amount of samples is present, this can be demanding (as the machine needs to keep a great number of very large files open at the same time)
  • Combining several mono-sample gVCF files into fewer multi-sample files can greatly speed-up the genotyping and be less resource demanding.

gatk CombineGVCFs \
	-O family.gvcf.gz \
	-R ref.fa \
	-V child.gvcf.gz \
	-V mother.gvcf.gz \
	-V father.gvcf.gz 

bash

3.5. Genotyping gVCF files

  • This step looks at the genotype of each sample at each position
  • It produces biallelic or multiallelic variants and affects genotypes to each sample
  • Only variant sites are kept
  • Several information (statistics, local frequencies, quality, etc.) are also provided.
gatk GenotypeGVCFs \
	-O variants.vcf.gz \
    -R ref.fa \
    --dbsnp /DATA/bioinfo/chr12M.dbSNP.146.b37.vcf.gz \
    -V family.gvcf.gz

bash

3.6. Recalibrating Variant Quality Scores

Calculates a new quality score called VQSLOD (Variant Quality Score Log-Odds)

  • Allowing to finely balance sensitivity (less FNs) and specificity(less FPs)
  • Uses Machine Learning to affect a score based on various contextual information
  • Trains on highly validated variant data (omni, 1000g, hapmap)
  • ⚠ Large amount of data is needed to perform this step! At least exome wide and with as many samples as possible
  • The user defines a set of quality tranches, containing a certain amount of True/False Positives
  • Later, these tranches will be used to filter out data, depending on the needs.

3.7. VCF File Format

VCF

A VCF File, is a text file consisting of a Header (starting with #) and record lines. Each line refers to a genomic position bearing one or more variants.

There should only be a line if at least one sample carries one of the variant alleles.

The data in the lines are distributed over several columns :

# Column Description Examples
1 CHROM the chromosome of the variant site 12, chr17, X, chrY
2 POS the start position of the variation 156701805
3 ID an identifier for the variant empty, ., rs179811, anything
4 REF the allele in the Reference Genome at the given position A
5 ALT the list of alternate alleles present in the samples C,G
6 QUAL the variant calling quality 1256.3
7 FILTER a list of reasons why we should ignore this variant ., PASS, VQSTRTranche90to99
8 INFO Metrics and annotation about this variants
9 FORMAT List of fields contained in the following samples columns GT:DP:AD:PL:GQ
10+ samples one column per sample, containing the genotype and other values 0/1:33:16,17:65,0,870:65
Allele Numbering and Genotype Notation

Allele Numbers

The Reference allele REF is numbered 0
The Alternate alleles ALT, ALT,...,ALT are numbered 1, 2, n

Example with REF=A and ALT=C,G

  • 0 ➡A
  • 1 ➡C
  • 2 ➡G
  • . ➡ unknown

Genotypes

Human having 2 chromosomes in each pair, 2 values are expected. Those values are separated by /or |

Example with REF=A and ALT=C,G

  • 0/0 ➡ A/A
  • 0/1 ➡ A/C
  • 1/1 ➡ C/C
  • 1/2 ➡ C/G
  • ./. ou . ➡ missing
Phasing

/ indicates unphased genotypes and | indicates phased genotypes

  • if a given individual has a variant A/C, then a variant G/T ➡this individual has an allele A and C, then an allele G then T
  • if a given individual has a variant A|C, then a variant G|T ➡ alleles A and G are on the same chromosome and alleles C and T are on the other chromosome. Thus AG and CT form haplotypes
How to

To answer the following question we'll examine the variants.vcf.gz file. We'll use the commands seen previously (zcat, cut, grep, wc, ...)

3.7.1. What is the version of this VCF file ?

Look at the fileformat header.

use grep

3.7.2. Which version of GATK was used to produce this files ?

Look for GATKCommandLine in the header.

use grep

3.7.3. What is the meaning of the following keyword ?

  • AC (defined in the INFO header)
  • GT (defined in the FORMAT header)
  • DP (defined in the FORMAT header and the INFO header)
  • GQ (defined in the FORMAT header and the INFO header)

use grep

3.7.4. List the individuals (samples) in this file

use grep, tail and cut

3.7.5. How many variants are there in this file ?

That is:

  • Number of variants sites ➡easy
  • Number of variants alleles ➡ harder

use grep -v, wc (and cut, tr)

3.7.6. How to differentiate SNVs, insertions, deletions and multiallelic variants ?

Look at the variants starting with C:

zcat variants.vcf.gz | grep -v "^#" | cut -f4,5 | grep "^C" | sort -u | column -t

bash

3.7.7. Examine the following variants

  • 12:11214360 (rs73260771)
    • How many samples have the following alleles: A, C and T ?
    • What is the allele count for the alleles A, C and T ?
  • 12:9833524 (rs3764022)
    • How many samples are homozygous / heterozygous ?
    • What is the GQ values for each sample ?
    • What is the DP values for each sample ?
  • 12:129298819 (rs3741614)
    • How many alternates are substitution, insertion, deletion ?
    • How many genotypes are phased / unphased ?

Next Step: 4 Databases and Viewers